home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / program / nrpas13.zip / EL2.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  1KB  |  57 lines
  1. FUNCTION el2(x,qqc,aa,bb: real): real;
  2. LABEL 1;
  3. CONST
  4.    pi=3.14159265;
  5.    ca=0.0003;
  6.    cb=1.0e-9;
  7. VAR
  8.    a,b,c,d,e,f,g: real;
  9.    em,eye,p,qc,y,z: real;
  10.    l: integer;
  11. BEGIN
  12.    IF (x = 0.0) THEN el2 := 0.0
  13.    ELSE IF (qqc <> 0.0) THEN BEGIN
  14.       qc := qqc;
  15.       a := aa;
  16.       b := bb;
  17.       c := sqr(x);
  18.       d := 1.0+c;
  19.       p := sqrt((1.0+c*sqr(qc))/d);
  20.       d := x/d;
  21.       c := d/(2.0*p);
  22.       z := a-b;
  23.       eye := a;
  24.       a := 0.5*(b+a);
  25.       y := abs(1.0/x);
  26.       f := 0.0;
  27.       l := 0;
  28.       em := 1.0;
  29.       qc := abs(qc);
  30. 1:      b := eye*qc+b;
  31.       e := em*qc;
  32.       g := e/p;
  33.       d := f*g+d;
  34.       f := c;
  35.       eye := a;
  36.       p := g+p;
  37.       c := 0.5*(d/p+c);
  38.       g := em;
  39.       em := qc+em;
  40.       a := 0.5*(b/em+a);
  41.       y := -e/y+y;
  42.       IF (y = 0.0) THEN y := sqrt(e)*cb;
  43.       IF (abs(g-qc) > ca*g) THEN BEGIN
  44.          qc := sqrt(e)*2.0;
  45.          l := l+l;
  46.          IF (y < 0.0) THEN l := l+1;
  47.          GOTO 1
  48.       END;
  49.       IF (y < 0.0) THEN l := l+1;
  50.       e := (arctan(em/y)+pi*l)*a/em;
  51.       IF (x < 0.0) THEN e := -e;
  52.       el2 := e+c*z
  53.    END ELSE BEGIN
  54.       writeln('pause in routine EL2'); readln
  55.    END
  56. END;
  57.